# Author: Mark Richardson
# Purpose: Estimate alternative model specifications and conduct sensitivity analysis
# Appendix F

# Load packages
library(dplyr)
library(fixest)
library(sensemakr)
library(kableExtra)
library(pander)

# Load data
perf <- readxl::read_excel(path = "data/01_lpr_012824.xlsx")

# create variables
perf <- perf %>%
  mutate(monthsVacant = (1 + Total_Vacant_Days) / 30,
         ObamaMvac = (1 + Obama_Vacancy_Length_Total) / 30,
         BushMvac = (1 + Bush_Vacancy_Length_Total) / 30,
         emp = employ16 / 1000)

# Rename variabels for sensemakr table
perf <- perf %>%
  rename(agencyPerformance = hier_perf_in,
         workforceSkill = skills_mean_2014,
         agencyIdeology = ideo_rating)

#### Table F1 ####

m1 <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
              cabinet + offsec + indcom + depprior + workforceSkill +
              agencyIdeology,
            data = perf,
            vcov = ~dept_acr)

m1_abs <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
                  cabinet + offsec + indcom + depprior + workforceSkill +
                  abs(agencyIdeology),
                data = perf,
                vcov = ~dept_acr)

m1_sq <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
                 cabinet + offsec + indcom + depprior + workforceSkill +
                 agencyIdeology + I(agencyIdeology^2),
               data = perf,
               vcov = ~dept_acr)

m1_lag <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
                  cabinet + offsec + indcom + depprior + workforceSkill +
                  abs(agencyIdeology) + ObamaMvac,
                data = perf,
                vcov = ~dept_acr)

m1_lag_2 <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
                    cabinet + offsec + indcom + depprior + workforceSkill +
                    abs(agencyIdeology) + BushMvac,
                  data = perf,
                  vcov = ~dept_acr)



setFixest_dict(agencyPerformance = "Agency Performance Rating",
               monthsVacant = "Months Vacant",
               directPAS = "Direct PAS Leader",
               totturn = "Leadership Transitions",
               eop = "EOP",
               cabinet = "Cabinet Department",
               offsec = "Whole Department",
               indcom = "Independent Commission",
               depprior = "Priority Department",
               workforceSkill = "Workforce Skill -- Obama Admin.",
               agencyIdeology = "Agency Ideology",
               `abs(agencyIdeology)` = "|Agency Ideology|",
               `I(agencyIdeology^2)` = "Agency Ideology$^2$",
               ObamaMvac = "Total Mths Vacant -- Obama Admin.",
               BushMvac = "Total Mths. Vacant -- Bush (43) Admin.",
               dept_acr = "Dept.",
               strg_pty_dis = "Intense Partisan Disagreement",
               is_pres_priority = "Presidential Priority")


# Table in the appendix was created in Rmarkdown with the tex argument set to
# TRUE
etable(m1, m1_sq, m1_abs, m1_lag, m1_lag_2,
       tex = FALSE,
       digits = "r3",
       digits.stats = "r2",
       order = "!Constant",
       style.tex = style.tex(stats.title = "\\midrule"),
       title = "Alternative Model Specifications",
       label = "table1")

## 
## # Value of ideology that maximizes performance
## -0.051548  / (2 * -0.087046)
## 

#### Table F2 ####

## Data not provided for models in Table F2

## m_pr <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
##               cabinet + offsec + indcom + workforceSkill +
##               agencyIdeology + is_pres_priority,
##            data = perf,
##            vcov = ~dept_acr)
## 
## m_pr_null <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
##               cabinet + offsec + indcom + workforceSkill +
##               agencyIdeology,
##            data = perf %>% filter(!is.na(is_pres_priority)),
##            vcov = ~dept_acr)
## 
## m_dis <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
##               cabinet + offsec + indcom + depprior + workforceSkill +
##               strg_pty_dis,
##            data = perf,
##            vcov = ~dept_acr)
## 
## m_dis_null <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
##               cabinet + offsec + indcom + depprior + workforceSkill,
##            data = perf %>% filter(!is.na(strg_pty_dis)),
##            vcov = ~dept_acr)
## 
## m_bth <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
##               cabinet + offsec + indcom + is_pres_priority + workforceSkill +
##               strg_pty_dis,
##            data = perf,
##            vcov = ~dept_acr)
## 
## m_bth_null <- feols(agencyPerformance ~ monthsVacant + directPAS + totturn + eop +
##               cabinet + offsec + indcom + is_pres_priority+ workforceSkill,
##            data = perf %>% filter(!is.na(strg_pty_dis) & !is.na(is_pres_priority)),
##            vcov = ~dept_acr)
## 
## 
## 


## Table F2
## Table in the appendix was created in Rmarkdown with the tex argument set to
## TRUE
## etable(m_pr, m_pr_null, m_dis, m_dis_null, m_bth, m_bth_null,
##        tex = TRUE,
##        digits = "r3",
##        digits.stats = "r2",
##        order = "!Constant",
##        style.tex = style.tex(stats.title = "\\midrule"),
##        title = "Alternative Measures of Presidential Priority and Partisan Disagreement",
##        label = "table2")
##

## Table F3 - see 01_r_p_l_jop_012724_del.do

#### Sensitivity Analysis ####

# Note: Homoskedastic SE is larger than clustered SE

perf <- perf %>%
  rename(deptPriority = depprior)

m1 <- lm(agencyPerformance ~ monthsVacant + directPAS + totturn + eop + cabinet
         + offsec + indcom + deptPriority+ workforceSkill + agencyIdeology,
         data = perf)

m1_mv <- lm(agencyPerformance ~ directPAS + totturn + eop + cabinet + offsec
            + indcom + deptPriority + workforceSkill,
            data = perf)

# Main specification - all data
d1 <- lm(monthsVacant ~ directPAS + totturn + eop + cabinet + offsec + indcom +
           deptPriority + workforceSkill + agencyIdeology,
         data = perf)


m1_sense_skl <- sensemakr(model = m1,
                          treatment = "monthsVacant",
                          benchmark_covariates = "workforceSkill",
                          kd = c(1),
                          ky = c(1), # Strength of confounder in terms of benchmark
                          q = 1,
                          alpha = 0.05,
                          reduce = TRUE)

m1_sense_pr <- sensemakr(model = m1,
                         treatment = "monthsVacant",
                         benchmark_covariates = "deptPriority",
                         kd = c(1),
                         ky = c(1), # Strength of confounder in terms of benchmark
                         q = 1,
                         alpha = 0.05,
                         reduce = TRUE)

m3 <- lm(agencyPerformance ~ monthsVacant + directPAS + totturn + eop + cabinet
         + offsec + indcom + deptPriority+ workforceSkill + abs(agencyIdeology),
         data = perf)

m3_sense_id <- sensemakr(model = m3,
                         treatment = "monthsVacant",
                         benchmark_covariates = "abs(agencyIdeology)",
                         kd = c(1),
                         ky = c(1), # Strength of confounder in terms of benchmark
                         q = 1,
                         alpha = 0.05,
                         reduce = TRUE)

m4 <- lm(agencyPerformance ~ monthsVacant + directPAS + totturn + eop + cabinet
         + offsec + indcom + deptPriority+ workforceSkill + abs(agencyIdeology) + ObamaMvac,
         data = perf)

m4_sense_ob <- sensemakr(model = m4,
                         treatment = "monthsVacant",
                         benchmark_covariates = "ObamaMvac",
                         kd = c(1),
                         ky = c(1), # Strength of confounder in terms of benchmark
                         q = 1,
                         alpha = 0.05,
                         reduce = TRUE)

# Table F4 - see the list of three items under "Sensitivity Statistics"

m1_sense_skl 

# Table F5

# Make a tibble to pass to kbl()

bench <- tibble(var = c("Workforce Skills", "Priority Department", "|Agency Ideology|", "Months Vacant (Obama)"),
                r2yz_dx = c(m1_sense_skl$bounds$r2yz.dx,
                            m1_sense_pr$bounds$r2yz.dx,
                            m3_sense_id$bounds$r2yz.dx,
                            m4_sense_ob$bounds$r2yz.dx) * 100,
                r2dz_x = c(m1_sense_skl$bounds$r2dz.x,
                           m1_sense_pr$bounds$r2dz.x,
                           m3_sense_id$bounds$r2dz.x,
                           m4_sense_ob$bounds$r2dz.x) * 100,
                r_qa = c(m1_sense_skl$sensitivity_stats$rv_qa,
                         m1_sense_pr$sensitivity_stats$rv_qa,
                         m3_sense_id$sensitivity_stats$rv_qa,
                         m4_sense_ob$sensitivity_stats$rv_qa) * 100,
                r2y_dx = c(m1_sense_skl$sensitivity_stats$r2yd.x,
                           m1_sense_pr$sensitivity_stats$r2yd.x,
                           m3_sense_id$sensitivity_stats$r2yd.x,
                           m4_sense_ob$sensitivity_stats$r2yd.x) * 100)

# Create the table (formatted for Rmarkdown)

## kbl(bench, booktabs = TRUE, digits = 2,
##     col.names = c("Benchmark Variable",
##                   "$R^2_{Y \\sim Z | \\textbf{X},D}$",
##                   "$R^2_{D \\sim Z | \\textbf{X}}$",
##                   "$RV_{q=1,\\alpha = 0.5}$",
##                   "$R^2_{Y \\sim D | \\textbf{X}}$"),
##     escape = FALSE, label = "benchTble",
##     caption = "Benchmark Confounders",
##     align = "lcc|cc") %>%
##   kable_styling(latex_options = "striped")
## 

kbl(bench, digits = 2)
